Majorana state on the surface of a disordered 3D topological insulator 
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We study low-lying electron levels in an "antidot" capturing a coreless vortex on the surface of 
a three-dimensional topological insulator in the presence of disorder. The surface is covered with 
a superconductor film with a hole of size R larger than coherence length, which induces supercon- 
ductivity via proximity effect. Spectrum of electron states inside the hole is sensitive to disorder, 
however, topological properties of the system give rise to a robust Majorana bound state at zero 
energy. We calculate the subgap density of states with both energy and spatial resolution using the 
supersymmetric sigma model method. Tunneling into the hole region is sensitive to the Majorana 
level and exhibits resonant Andreev reflection at zero energy. 
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Topological insulators and superconductors are very 
peculiar materials with a gap in the bulk electron spec- 
trum and a low-lying branch of subgap excitations on 
their surface (see [l| for a review). This surface metallic 
state appears due to topological reasons and is robust 
with respect to any (sufficiently small) perturbations. In 
particular, topological properties prevent these surface 
states from Anderson localization. One common exam- 
ple of a topological insulator is a two-dimensional system 
in the integer quantum Hall effect regime. The bulk of 
such a system has a spectral gap between successive Lan- 
dau levels and is hence an insulator. At the same time 
quantized Hall conductance appears due to a fixed in- 
teger number of chiral propagating edge modes on the 
background of the bulk gap. This type of materials are 
referred to as Z topological insulators. 

Another type of topological insulator is realized in 
three-dimensional (3D) semiconductor crystals with suf- 
ficiently strong spin-orbit interaction (BiSb, BiTe, BiSe, 
strained HgTe etc). The spin-orbit interaction leads to 
inversion of the spectral gap. As a result subgap surface 
excitations appear with a dispersion of the massless Dirac 
type. The topological invariant in these materials has a 
Z2 nature. When the number of surface states is odd, 
one of them always remains gapless due to topological 
protection. 

A general classification of topological insulators was 
developed in Refs. 0, 0| based on the symmetry of the 
underlying Hamiltonian. The quantum Hall effect is an 
example of a topologically nontrivial state of the uni- 
tary symmetry (class A of the Altland-Zirnbauer clas- 
sification 0|) in 2D. Strong spin-orbit interaction leads 
to a topological state in the symplectic class (All) in 
3D. Another important example is the ID topological 
superconductor of the class BD symmetry (superconduc- 
tor with both time-reversal and spin rotation symmetries 
broken). The topologically protected mode in this case 
is zero-dimensional and is known as the Majorana bound 
state (MBS). It appears, in particular, in the core of an 
Abrikosov vortex in the spinless p-wave superconductor 
A similar MBS appears 0] in a vortex in an ordinary 
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Figure f : (Color online) A sketch of the considered setup. 
Three-dimensional topological insulator is covered by a su- 
perconducting film with a hole of radius R. The gap in the 
density of states is induced on the covered surface forming 
an antidot and confining surface excitations inside the hole. 
External magnetic field induces an Abrikosov vortex inside 
the hole. The spectrum of energy eigenstates inside the hole 
acquires a Majorana zero level which can be accessed in the 
tunneling measurement. 



s-wave superconductor brought in contact with the sur- 
face of a 3D Z2 topological insulator, which corresponds 
to the symmetry class AIL The MBS appears as a descen- 
dant of the topologically protected massless Dirac state 
on the free surface of the topological insulator. 

There are two general methods to observe the Majo- 
rana level. One of them relates to the anomalous Joseph- 
son effect Another, and a more direct, way involves 
tunneling into the region where the MBS is supposed to 
be localized [9]. The differential conductance in such a 
tunneling experiment yields the local density of states 
with spatial and energy resolution. The Majorana state 
in the vortex core manifests itself by resonant Andreev 
reflection at zero energy [Icj . 

In this Letter we study local tunneling conductance 
for a setup depicted in Fig. [TJ A superconducting (with 
a spectral gap A) film with a circular hole of the radius 
R is deposited on the surface of a 3D topological insu- 
lator, e.g., Bi2Te3 crystal. Perpendicular magnetic field 
applied to the system produces a vortex pinned to the 
hole. The radius R is supposed to be relatively large, 
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R > (l,£o), where I is the mean free path for topologi- 
cal insulator surface states and £o is the superconducting 
coherence length of the film. The condition R > £ al- 
lows one to avoid the abundance of low-lying Caroli- 
de Gennes-Matricon states whose presence in the core 
of the Abrikosov vortex complicates observation of the 
E = Majorana state. Under these conditions disorder 
is relevant for the states localized in the hole region. Al- 
though the Majorana state is protected against disorder, 
i.e. its energy stays exactly zero, the spatial distribution 
of its wavefunction ipo(r), and thus the local tunneling 
conductance are sensitive to disorder. The aim of this 
Letter is to calculate the tunneling conductance in the 
"dirty-limit" with spatial and energy resolution, and to 
identify the effects of the Majorana zero-energy state in 
the presence of strong disorder. 

The problem is described by the following Bogolyubov- 
de Gennes (BdG) Hamiltonian in polar coordinates r, ip: 



H = 



A(r)e i 
A(r)e"^ -H 



H = v Q s • p + V(r) - fi. 

(1) 

Here the Hamiltonian Hq describes the dynamics of sur- 
face excitations in a topological insulator without a su- 
perconducting layer. The Fermi velocity of surface elec- 
trons is denoted by vq, s is the spin operator, and V(r) is 
the random disorder potential. The Fermi level is shifted 
from the Dirac point by fi. The vector potential term in 
Hq is neglected due to smallness of magnetic field. We 
assume the dirty limit with a disorder induced mean free 
path I <C £ = y/HD/A and a relatively large hole with 
R ^> £; here D = Vqt/2. These conditions allow us to 
use a step-like radial dependence of the order parameter: 



A(r) 




(2) 



Inside the hole, the order parameter is zero and elec- 
tron dynamics is governed by Hq. This Hamiltonian 
possesses time-reversal symmetry of the symplectic type, 
Hq = SyHqSy, and hence belongs to the symplectic sym- 
metry class AIL Proximity to the superconductor induces 
a gap in the electronic spectrum outside the hole. This 
gap effectively confines the low-lying excitations and im- 
poses boundary conditions for the Hamiltonian Hq at 
r = R. These boundary conditions break time-reversal 
symmetry due to the spatially rotating phase of the order 
parameter. At the same time the total BdG Hamiltonian 
H acquires a specific particle-hole symmetry: 



H 



-SyTyH SyTy. 



(3) 



Here r y is the Pauli matrix acting in Nambu-Gor'kov 
(NG) space. Identity ((3]) implies a symmetry of the elec- 
tron spectrum. For each eigenstate ip with energy E there 
is a conjugate eigenstate s v T y ip* with energy — E. 



The particle- hole symmetry Q defines the class BD. 
On the level of random matrices there are two distinct 
versions of this class with even and odd number of eigen- 
states referred to as D and B class, respectively. In class 
B, one unpaired eigenstate has exactly zero energy and 
is self-conjugate: ip — s y T y ip* . The BdG Hamiltonian 
counts every physical excitation twice due to the dou- 
bling in NG space. An unpaired level is thus "half" of a 
true excitation — a Majorana state. 

We will calculate the density of states inside the hole 
with the help of the supersymmetric non-linear sigma 
model. Two-dimensional Dirac fermions with potential 
disorder are described by a very peculiar model of the 
class All with Z2 topological term 



12]. This topologi- 



cal term appears as a consequence of the chiral anomaly 
of Dirac fermions. We will consider the minimal model 
operating with the 8x8 supermatrix Q. Apart from 
Fermi-Bose superspace, this matrix operates in the space 
of retarded-advanced (RA) fields and in a specific time- 
reversal (TR) space introduced to take into account 
Cooperons and diffusons on equal footing. RA space is 
completely analogous to the NG space (see below); we 
denote Pauli matrices in this space by r. Notation a is 
used for Pauli matrices in TR space. A detailed deriva- 
tion of the sigma model can be found in jl3| . 

In class All, the matrix Q obeys the non-linear con- 
straint Q 2 = 1 and the linear constraint 



Tr<T 



Q = Q = CQ 1 C 



C 



7 X 

icr. 



(4) 
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As a result, Q contains 8 commuting and 8 anticommut- 
ing (Grassmann) variables in total. Commuting degrees 
of freedom parameterize diagonal fermion-fermion (F) 
and boson-boson (B) blocks of Q. The general Q-matrix 
can be decomposed as 



Q = U- 



2f 
Qi 



U, 



(5) 



where the central part contains only commuting variables 
while Grassmann parameters define the unitary super- 
matrix U. We will use the following explicit form of the 
central part of Q in terms of eight angle parameters: 

Qf = [i~ z cosOf + a z sm9f(T x cos(f>f + r y sin </>/)] 

x [a z coskf + T z smkf(a x cosxf + cr y sinxf)], (6) 

Qb = r z cos 9b[cr z cos fa + sin h(a x cos Xb + o y sin Xb)] 
+ sin 9b(T x cos tfib + T y sin (pi,). (7) 

This representation fulfills all the constraints imposed by 
the symmetry class All of the Hamiltonian Hq. The F 
and B sectors of the sigma model are compact and non- 
compact, respectively. This is achieved by demanding 
that the angles 9b, kb are imaginary while all other angles 
in Eqs. ©, © are real. Below we will find that the saddle 
point describing the density of states in the hole occurs 
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with real 6b- This implies a proper shift of the integration 
contour for this angle. 

The sigma model of class All is designed to study 
transport properties of a disordered system. This im- 
plies that the Q matrix operates, in particular, in the 
RA space allowing for averaging the product of retarded 
and advanced Green functions. We are interested in the 
density of states and hence it suffices to average just the 
single retarded function. At the same time, supercon- 
ducting boundary conditions implemented in the BdG 
Hamiltonian ([1]) require to introduce an additional dou- 
bling of fields in the Nambu-Gor'kov space. This can be 
achieved within the standard All class sigma model with 
the 8x8 supermatrix Q while the role of NG space is 
taken over by the RA structure of Q (for detailed discus- 
sion of the transformation from RA to NG representation 

HQ). 



see 



Thus we can incorporate the superconducting order 
parameter directly into the action of the sigma model, 

S[Q] = ^ Jd 2 r Str [D{VQf + A(ieA - A)q] + S e [Q] , 

A = r z a z , A = A(r)(r x costp - t v simp). (8) 

Here e = E + iGtS(r — t^/A-kv is the sum of the energy 
E and the local dwell term describing the coupling to a 
tunnel tip with dimensionless conductance Gt <C 1 fl3j |. 

The action (JSJ involves the Z 2 topological term Sg[Q] 
which appears due to massless Dirac nature of underly- 
ing electrons as was discussed above. The topological 
term involves only the compact part of the sigma-modcl 
manifold, i.e., its F sector VWf- In the general version 
of the sigma model, that is capable of averaging several 
retarded and advanced Green functions, the homotopy 
group is 7T2(.Mf) = ■ However, in the minimal model 
we are considering, A^f has the structure of the product 
of two spheres S 2 x S 2 as seen from Eq. ([6|). In this case 
the homotopy group is richer, ZxZ. Two integer topolog- 
ical invariants can be introduced counting the degree of 
covering of the two spheres by the mapping Q from real 
space to A4p. This allows us to write the topological 
term explicitly although in a non-invariant form. With 
the parameterization ((6]), it reads (cf. Ref. [lH): 



Se[Q] = \j* 



sin0/(V0/ x V(/>/) 

+ sin£;/(Vfc/ x V\f) 



(9) 



Let us now analyze the minima of the action (JSj) for 
the setup Fig. [T] Circular symmetry of the problem fixes 
the phase equal to the polar angle, (j)f.b = ip. The other 
parameters depend only on the radial coordinate r. Both 
angles 8 in F and B sectors obey the Usadel equation 



D 



dr 2 



1 86 

r dr 



sin 26 



2r 2 

2iEsin(9cos£; + 2A(r)cos6l = 0. (10) 



Inside the hole A = 0. At low energies we can also neglect 
the E term and the equation becomes independent of k. 
The step- like dependence of the order parameter, Eq. @, 
imposes the boundary condition 6(R) = ir/2. There are 
two possible solutions to the Usadel equation with this 
boundary condition: 



91 = 2arctan(r/i?), 

9 2 = 7r — 2 arctan(r/i?). 



(11) 
(12) 



Saddle point equations also require kf — and hence the 
angle \f drops from the matrix Q and from the action. 
The two remaining angles kf, and Xb are free and can take 
any constant values. 

The spatial profile of Q is thus fixed by the Usadel 
equation. The solutions 0i ; 2 represent two disconnected 
saddle points in the F sector while in the B sector only 
the saddle point 6b — 9\ is reachable. If the integration 
contour for 6b, which runs along the imaginary axis, is 
shifted to the point 62 a divergence occurs in the kb in- 
tegral. Thus the B sector is reduced to a hyperboloid 
parameterized by ikb > and < Xb < 27r, see Ref. fDij . 
This is exactly the structure of the sigma model of class 
BD as we anticipated from symmetry analysis. The dis- 
tinction between even (D) and odd (B) versions of this 
class is related to the disjoint character of the manifold 
due to the discrete degree of freedom in the F sector. 
Namely, in class D (B) the two parts of the manifold 
contribute to the partition function with the same (op- 
posite) sign. In our problem the odd symmetry class B 
occurs. To demonstrate it, we compare the value of the 
action (j8]) at the two minima in the F sector. These two 
solutions indeed contribute with the opposite sign since 
the corresponding values of the topological term ([9]) differ 
by exactly iir. 

The density of states is given by the integral 



p(E,r) 



Re / DQStr[fcAQ(r)] e 



-S[Q] 



(13) 



At low energies, this integral is to be calculated over the 
saddle manifold. Apart from the two variables kb and Xb 
and two disconnected points in the F sector this manifold 
involves two Grassmann variables in the matrix U, Eq. 
(|5|). In order not to spoil the saddle point, these vari- 
ables must be constant in space and fulfill the condition 
[U, A] = 0. We introduce them according to 



U = exp 



Wx + (<7y 



FB 



(14) 



Within this parameterization, we can rewrite the integral 
(|13p in terms of kb, Xb, f]i C- Thus the problem is reduced 
to the 0D sigma model of class B (the two disjoint points 
in the F sector contribute with opposite signs). Explicit 
calculation of the integral [ljl, [HI yields the local DOS 
in the factorized form 



p{r,E) = vn{r)f{E/u ), 



(15) 
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Figure 2: (Color online) Global density of states as a func- 
tion of energy E. The solid blue oscillating curve represents 
the low-energy result (|19[l . the dashed curve shows the high- 
energy asymptotics (|20|) , and the red curve presents numerical 
solution interpolating between the two limits. 



/(*) = 



n(r) = cos 9\ (r) 

7 



R 2 ~r 2 
R 2 + r 2 ' 
sin(27rx) 



tt(x 2 + 7 2 ) 



1 - 



2nx 



(16) 
(17) 



Here 7 = G t n{ro)/2i: <C 1 and the low-energy level spac- 
ing is given by 

cj- 1 =2i> I d 2 r cos 0i (r) = 27r(log 4 - 1)jaR 2 . (18) 



The spatial profile of DOS, n(r), is fixed by the solution of 
the Usadel equation, while the energy dependence is char- 
acteristic for the B class and contains a narrow lorentzian 
peak at zero energy. This peak is the zero-energy MBS 
broadened due to the finite tunneling time. Note that 
the width 7 is position-dependent. 

Integrating p(r, E) over space yields the global DoS 



N(E) 



f(E/uo) 
2uj 



(19) 



In the extreme tunneling limit 7 — > 0, this function ac- 
quires the contribution S(E)/2 and coincides with the 
result [15[ up to a factor 2 due to BdG double counting. 
The Majorana state appears as a half of a fermionic level. 

Spatial and energy dependence of p(r, E) factorize at 
energies much less than the Thouless energy E^ = 
D/R 2 . At higher energies fluctuations of Q are not im- 
portant and DOS is given just by a single saddle point. 
Approximate solution of the Usadel equation (fT0|) in the 
limit E > E Th yields [H 



N{E > E Th ) = tivR 2 1 - (2 - V2) ^E Th /E . (20) 

Local DOS is close to the normal value v everywhere 
except for a narrow vicinity of the hole boundary. 



The local density is measured in the tunneling experi- 
ment. The tunneling current is determined by 

I= ^tj P( E ' r °) U( E ~ eV ) ~ dE (21) 

with f(E) being the equilibrium Fermi distribution func- 
tion. At low temperatures and voltages, (T, eV) <C ujq, 
we keep only the first lorentzian term in the energy de- 
pendence of the local DOS (p~7|) . The differential conduc- 
tance exhibits a peak 



eV 



dl_ = J 7rfi.[ 7 2 + (eV7u;o) 2 ] 
dV 

' AKT cosh 2 (eV/2T)' 



T < 7CJ0, 



7^0 « T < u . 



(22) 



At very low temperatures the height of this peak is uni- 
versal and equals e 2 /irh. This signifies resonant Andreev 
reflection at the Majorana state. This effect in the ab- 
sence of disorder was studied in Ref. [10] . At larger (and 
more realistic) temperatures, 7W0 <C T < cjo, the height 
of the peak is parametrically small, dl/dV ~ jljo/T. 

Noise power of the tunneling current in the same 
regime (T, eV) < w is [HI 



S(V,T,r ) 



e 2 ju>o 
2h : 



-fujQ <C T <C wo- 



(23) 



The noise produced by the Majorana level is T- and V- 
independent as long as jujq <C T. 

The results (|2^|) and (|23p apply at low temperature and 
voltage. When temperature and/or voltage are higher 
than the level spacing uiq, non-zero-energy states also 
contribute to the tunneling current. Positions and widths 
of these states depend on the realization of disorder. For 
low temperature but high voltage, T < ujq < eV, narrow 
resonances similar to Eq. (|22p will occur due to non-zero 
levels [l|| . Positions of these resonances strongly depend 
on disorder realization; their heights are smaller but close 
to e 2 /ttH and widths are of order 7^0- When tempera- 
ture exceeds uio, all the resonances get smeared and the 
normal average DOS v is recovered. 

To conclude, we have studied the local density of states 
in the superconducting vortex on the surface of a topo- 
logical insulator in the superconducting antidot setup de- 
picted in Fig. [TJ The spatial and energy dependence of 
the density of states factorize at low energies and the lat- 
ter is given by the 0D sigma model of symmetry class B. 
We have identified the zero-energy Majorana state occur- 
ring due to the topological properties of the system. This 
Majorana state exhibits itself via the resonant Andreev 
reflection at zero energy yielding the peak in differen- 
tial conductance with the universal amplitude e 2 /irH and 
width proportional to the normal conductance Gt- 
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Supplementary Information 
Derivation of the sigma-model 

Sigma model for hybrid structure 



In this section we derive the sigma model describing the dynamics of the topological insulator surface excitations in 
the presence of a random potential. We assume the hole geometry depicted in Fig.[T]of the main text. The derivation 
starts with the microscopic Bogolyubov - de Gennes Hamiltonian ([I]). The sigma model is aimed at averaging the 
single retarded Green function determining the density of states. We will also demonstrate the equivalence of this 
sigma model to the model of the symplectic symmetry class All describing the transport properties of massless Dirac 
fermions. The latter model is based on the normal, rather than superconducting, Hamiltonian and is capable of 
calculating averaged products of retarded and advanced Green functions. It is known that the sigma model for Dirac 
fermions possesses the specific Z 2 topological term. The equivalence of the two models will prove the appearance of 
the topological term in the superconducting case considered in the present paper. 

Let us start with the retarded Green function in the system governed by the Hamiltonian (JXJ) . We will use the 
supersymmctric integral representation 



G|(r,r') 



■J str [fc$(r)$t(r') 



5' 



d 2 r&(E + i0 - H)$. 



(24) 



Here the vector fields <J> and & contain 4 commuting and 4 Grassmann parameters each. Apart from the spin 
and Bogolyubov - de Gennes structure of the Hamiltonian, we also introduce the superstructure in order to get rid 
of the normalizing denominator and facilitate further disorder averaging. The pre-exponential factor contains the 
supermatrix k = {1, 



17]: str A = A 



FF 



A 



BB- 



1}fb and supertrace is defined as in Rcf. 
The superconducting symmetry of the Hamiltonian H, Eq. ©, gives rise to specific soft modes - Cooperons - of 
the type (G E G^ E ). These modes are relevant for the average density of states since they are built out of retarded 
functions only. In order to include them into our effective theory, we transform the action by writing half of it in the 
form (|2"4")1 and another half in the time-reversed (transposed) form: 



s 



d 2 r 



& (E + lO - H)<f> - <p T k(E + i0 + SyTyHSyTy)^ 

cPr (<I> + t,, i$ T fcs„r, 



(E + i0)a z T z -H -t z A 



8yTy& 



(25) 



Here we have introduced the matrix a z operating in the space of time-reversed blocks - TR space. The matrix k 
appeared due to anticommutation of Grassmann variables. We will denote the doubled vectors as 



* = 



1 



1 



* = — = (<S> + t z , i^ T ks v T x ) 
v2 



V2 \SyTy<$>* 

They are no longer independent, as $ and $* were, but rather obey the linear constraint 

'o x 



y * ia. 



(26) 



(27) 



y/ fb 
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At this stage, we are ready to perform the disorder averaging. We adopt the standard Gaussian white noise disorder 
model characterized by the correlator 



(V(r)V(r>)) = 



6(r 



(28) 



The single-particle density of states is defined as v = |/x|/27TUq. Note that the right-hand side is twice bigger compared 
to the standard definition applied to normal metals. The reason for this extra factor 2 is the Dirac nature of electron 
spectrum and hence doubling of the phase space due to the two types of excitations (electrons and holes). 
Disorder averaging produces the quartic term in the action, 



d 2 r 



1-KVT 



(E + iO)A + (j, - v s • p - t z A 



(29) 



Here we use the notation A = a z r z . The interaction term is further decoupled with the help of the Hubbard- 
Stratonovich transformation by introducing an auxiliary matrix field Q. The construction of the nonlinear sigma 
model implies that the field Q contains all relevant slow modes of the disordered system - diffusons and Cooperons. 
There are in total three ways to decouple the four-fermion interaction [T3] - One of them involves a scalar field ~ ('i n I') 
analogous to the random potential V. This leads only to an irrelevant renormalization of the chemical potential. Two 
other decoupling schemes introduce a matrix field with the structure (\E"I r ) and (^ r 4 rT ) corresponding to diffusons and 
Cooperons, respectively. In order to include all relevant slow modes in the theory, we have to perform decoupling in 
both ways, which can be achieved with a single matrix Q. Details of this calculation can be found in Ref. 17]. The 
resulting action has the form: 



16r 



strQ 2 



EA + fi — vq s ■ p — t z A 



iQ 
2r 



(30) 



The vectors $ and *S> are related by Eq. (|27|) . This allows us to limit the matrix Q by the linear constraint Q = Q = 
CQ T C T . This relation keeps only those parameters in Q that couple to the product 'J"!'. Integrating out the field $ 
yields the action for the matrix Q: 



16t 



strQ 2 



1 



str In 



EA + n — v s ■ p — t z A 



2r 



(31) 



Derivation of the sigma model proceeds with the saddle-point analysis of the above action. In the dirty limit, the 
energy and A terms in the argument of the logarithm are relatively small compared to Q/t. With these small terms 
neglected, the action possesses the uniform saddle point Q — A. This saddle point corresponds to the self energy of 
the average electron Green function in the self-consistent Born approximation. Other saddle points can be achieved 
by rotations Q = T —1 AT, if the matrix T commutes with the spin operator s (E and A terms are still neglected). 
Rotations T define the target manifold of the non-linear sigma model. The effective action of the sigma model within 
this manifold is derived with the help of the gradient expansion of Eq. pip allowing for slow spatial variation of Q 
and perturbative expansion to the linear order in E and A. Since the Q matrix is trivial in the spin space, we can 
safely reduce its size to 8 x 8 keeping only Nambu, TR, and FB structure. Then the self-conjugacy relation acquires 
the form of Eq. (gj). 

The gradient expansion of Eq. pip is a highly nontrivial procedure in view of the chiral anomaly of the Dirac 
fermions. The momentum integrals arising after the expansion of the logarithm are divergent in the ultraviolet limit 
and require a proper regularization. The result of the expansion is independent of a particular regularization scheme 
provided the gauge invariance is preserved. The anomaly affects only the imaginary part of the action and leads to 
the appearance of the topological term [l2[ ■ At the same time, the real part can be obtained in a straightforward way 
since all the arising momentum integrals are convergent. The result reads 



ReS=- 



d r str 



D(yQ) 2 +4i(EA-r z A)Q 



(32) 



Here the diffusion coefficient is D — v 2 t/2 and the matrix Q is reduced to the 8x8 size. The action ([5]), used in the 
main text, differs from Eq. (|32p only by a it/2 rotation of the superconducting phase and by an imaginary contribution 
to the energy term. The latter corresponds to a finite dwell time of the electron due to the coupling to the tunneling 
microscope tip as we elaborate below. 

In order to explain the emergence of the topological term in the imaginary part of the action, we will prove the 
equivalence of the sigma model, derived here for the hybrid structure of Fig. [TJ to the symplectic sigma model for 
Dirac fermions derived in Ref. 1121 . 
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Sigma model for Dirac fermions 



The sigma model obtained by the gradient expansion of Eq. (|31j) belongs to the symplectic symmetry class AIL 
This is quite natural since the Hamiltonian possesses only time-reversal symmetry in the absence of the A term. Let 
us demonstrate the equivalence between our sigma model (j8]) describing the density of states for the Hamiltonian H, 
Eq. ([I]), and the standard symplectic class sigma model describing transport via the surface states of the topological 
insulator with Hamiltonian H - In the latter case, the sigma model yields averaged products of retarded and advanced 
Green functions with the energy difference u>. This requires introducing retarded-advanced (RA) structure of the 
action, 

S o = -iJ d 2 r $V 2 [(w/2 + iO)r z - H ] $. (33) 

Here the matrix t 2 operates in the RA space. Subscript is used to distinguish this model from the sigma model 
derived in the previous section and used in the main text. 

The time-reversal symmetry H — s y H^ s y is taken into account by further doubling of variables in the TR space, 



J d 2 r $ f (ui/2 + iO - r 2 i7 )$ - $ T fc(u;/2 + iO - T z s y H s y )<f>* 

= -\\& 2 r ($tr 2 , i^ks y r z ) [(w/2 + *0)A - H ] . (34) 

Note that the matrix Ao = r 2 multiplying the frequency term is different from its counterpart A from the previous 
section. The vectors ^ and ^> as well as the charge conjugation matrix Co are also defined differently, 

* = ^ (4*) , * = «**) T = ± (*+r., iVkw.) , Co = -is v r z (- .°J FB . (35) 

Already at this stage we can prove the equivalence of the two models by selecting a proper unitary rotation of the 
field vector This unitary rotation should obey the following properties: 

U T C U = C, U^AqU = A. (36) 

These relations bring the charge conjugation matrix Co and the matrix Ao to the representation used in the previous 
section for the sigma model of the hybrid system. There are many matrices U that fulfill identities (f3f))) . One possible 
choice is 

In the case of usual rather than Dirac fermions, an analogous equivalence between the sigma model for the density of 
states in a superconducting hybrid structure and the orthogonal class sigma model is discussed in Ref. [bl ] . 



Topological term 

So far we have discussed the derivation of the real part of the sigma- model action ((5J) and also proved the equivalence 
of this model to the symplectic class sigma model for Dirac fermions (up to the boundary conditions involving the 
term A). The latter model is known to possess the 1^ topological term as an imaginary part of its action, see Ref. 
12;]. This allows us to include the topological term in the action © and thus complete its derivation. 

The explicit form of the topological term is written in the main text, Eq. @, in a noninvariant form, using the 
parameterization © - ([7]). In this section we will discuss an indirect but explicitly invariant representation of this 
topological term. 

The target manifold of the symplectic class sigma model is fixed by the representation Q = T _1 AT with the 
constraint Q = Q. If we neglect the Grassmann parameters and consider the central part of Q, see Eq. ([5]), these 
conditions yield Qf € 0(4)/0(2) x 0(2) and Qb G Sp(2, 2)/Sp{2) x Sp{2). The topological term arises in the compact 
F sector of the model. An explicit expression for the topological term can be written within a construction similar to 
the Wess-Zumino-Witten term, cf. Ref. [l8j . 
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We first extend the target manifold by relaxing the condition Q 2 = 1. Let us introduce the matrix Q — TAT. The 
only restriction on this matrix is Q = Q. The F sector of the extended manifold is Qf S 0(4). It has trivial second 
homotopy group, 7r 2 [0(4)] = 0. This allows us to introduce the third, auxiliary, coordinate < t < 1, such that the 
real 2D space corresponds to t = 1, and continuously extend the matrix field according to 

<*,.)-{«"■ (3 8 , 

Such an extension assumes that the physical space has no boundary and can be viewed as a surface of a 3D ball with 
t being the radial coordinate. With this definition of Q, the Wess-Zumino-Witten term has the form 

dt J d 2 rstr[Q- 1 (V a Q)Q- 1 {S7 b Q)Q- 1 {S7 c Q)\. (39) 

Here the indices a, b, and c take three values corresponding to the three coordinates t, x, and y, and e a 6 C is the unit 
antisymmetric tensor. 

The integrand in this expression (|39p is a total derivative therefore the result of the integration depends only on 
the value of Q at the boundary of integration domain, i.e., on the physical field Q = Q\t=i- This can be checked by 
an explicit calculation of the variation of Sq: 

SS e [Q} = 1 -^ J dt J rf 2 rV Q str[Q- 1 ( 5QQ- 1 (V b Q)Q- 1 (V c Q)] = JL J d 2 r str [q-HQ[Q- 1 (?7 x Q), Q- l (V v Q)] }. 

(40) 

A topological term does not change under continuous variations of the field Q but takes different values in different 
(disconnected) topological sectors of the model. The general Wess-Zumino-Witten term has a nonzero variation and 
hence does not obey this property. However, the term ()39j) is constrained by Q 2 = 1. Under this condition, the 
variation (|40|) is identically zero. Hence the Wess-Zumino-Witten term (|39| indeed plays the role of a topological 
term. 

The value of the topological term is quantized, yielding the topological charge of the field configuration, only in a 
system without boundary. This is not the case for the geometry considered in the paper. The matrix Q is defined in 
the finite region inside the hole (see Fig. [T] of the main text) with the boundary conditions fixed by the A term in the 
Hamiltonian (JXJ) - The value of the topological term is not completely fixed in this case. In particular, the expression 
(|39p depends not only on physical values of Q but also on the way it is extended in the third dimension, Eq. (|38l) . 
However, this uncertainty leads only to an uncontrolled imaginary constant in the action. This constant is the same 
in both topological sectors of the model and hence does not alter any observable quantities. 

An alternative derivation of the sigma-model action, including the topological term in the form Eq. (|39j) is possible 
within the non-abelian bosonization formalism. Let us for a moment assume, that the Hamiltonian has an extra chiral 
symmetry, H = —s z Hs z . This situation corresponds to the symmetry class Dili and can be realized, e.g., at the 
Dirac point in the presence of a random velocity disorder. The sigma model of the class Dili has an extended target 
space, corresponding to the manifold of Q introduced above. This sigma model possesses the Wess-Zumino-Witten 
term (|39|) . see Ref. [191 ], and is the result of the non-abelian bosonization of the initial fermionic problem. Non-zero 
chemical potential drives the system away from the Dirac point and breaks the chiral symmetry. This reduces the 
model to the symplectic symmetry class and restricts the field by the condition Q 2 = 1. The resulting action has the 
topological term in the form of the restricted Wess-Zumino-Witten term (|39l) as discussed above. 



Tunneling coupling term 

In this section we discuss the appearance of the imaginary contribution to the energy due to the presence of the 
tunneling tip. In order to include the coupling to an external metallic probe, we extend the Hamiltonian as 

" " (" L) ■ (41) 

Diagonal blocks of this matrix are the Hamiltonian H of the topological insulator surface ([T]) and the Hamiltonian 
Hm describing electron states in the metallic tip. Off-diagonal elements t and t' are tunneling amplitudes between 
the two subsystems. We assume these tunneling amplitudes are local in space and couple the states at point r$ on 
the sample surface to the states at the tip of the metallic probe. 
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Deriving the sigma model as described in previous sections we arrive at the expression similar to Eq. 
an extended matrix structure: 



16t 



strQ 2 



16t m 



str Q 



1 



M 



str In 



ea - t z u 



iQ/2r 
iQm/2t m 



f3B but with 



(42) 



The quantities vm, t~m, and Qm refer to the metallic tip. In the absence of tunneling, t = 0, the two subsystems 
decouple and we obtain independent sigma models describing topological insulator and the tip. We will assume the 
density of states vm to be large, that will fix Qm = A. Tunneling coupling is assumed to be weak. This allows us 
to expand the logarithm in powers of t and v. The lowest non- vanishing contribution appears in the second order of 
the perturbation theory and yields the tunneling term in the sigma model, 



• str 



(EA - t z H + iQ/2r)~ 1 T z t(EA - t z H m + 



v v M str (QT z tQ M T z t^) . 



(43) 



In the last expression we have used the saddle-point relation between Q and the Green function (EA — T z H + iQ/2r)~ 1 
in coincident points. We introduce normal dimensionless (in units e 2 /h) tunnel conductance of the junction Gt — 
tt 2 uvm tr(ti ;') (note that v and vm are total rather than spinless densities of states at the Fermi energy). Assuming 



Q 



M 



A, we rewrite the action as 



St = --^strAQ. 



(44) 



Thus the tunneling term has the same structure as the energy term and produces an imaginary contribution to E, 
see Eq. © of the main text. 



Calculation of the density of states 

Low energies E <C Eth 

The local density of states on the surface of the topological insulator is given by Eq. (fT3")) of the main text. The 
integral over matrix Q is to be calculated over the saddle manifold fixed by the solution of the Usadel equation (|10p . 
There are two distinct solutions 0\^, [see Eqs. (fTTj) . (fl"2"|) of the main text] for the angle 0-p. In the F sector of the 
model, all other parameters are fixed and thus 6*1,2 represent two disjoint parts of the integration manifold. At the 
same time, only the solution Ob — 0\ is allowed in the B sector, and two other angles, fce and xb, are free forming a 
hyperboloid (the parameter fee must be imaginary in order to ensure convergence of the integral) . 

To calculate the density of states from Eq. f(T3|) . we have to express Str(feAQ(r)), S[Q] and DQ in terms of kb, 
Xb, f]-, and £ for the two disjoint parts of the manifold. However, a direct use of our parameterization leads to an 
uncertainty in the integral over the part corresponding to Of — 82- The integral over drjd£ is zero, while the integral 
over kb diverges. We resolve this uncertainty by a trick proposed in Ref. [IBj . The two disjoint parts of the integration 
manifold can be mapped onto each other by the similarity transformation 

Q^T-'QT, T= ^ J) (45) 

We use the parametrization ([6]), {7J in the domain with 9f = 6% and then apply the transformation (|45|) to cover the 
second topological sector with Of — O2. Owing to the relation (|45|) . integration measure is the same in both sectors. 
It is given by the superdeterminant of the metric tensor in the space of Q matrices, see Ref. [T3|- The integration 
measure reads 

1 K 

DQ = — tanh — dn dxb dr] d(, (46) 

Z7T 2 

with k = ikb (kb is imaginary to ensure the convergence of the integral). The integration runs over the hyperboloid 
k > 0, < Xb < 2tt. 

The action in the two topological sectors Of = 0\ ; 2 has the form 

/2 ( 2 K 2 K \* 7r / , 2 /f , K\ 11 
d r e cos 0i I sinh — — zryCcosh — J — = —Imx I sinh" — — i?7Ccosh — J — , (47) 

d 2 re cos X cosh 2 - - — = -2i7ricosh 2 - + — . (48) 
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Here we have introduced the complex dimensionless energy parameter x = J d 2 r e cos Q\ — x + iy, with the real part 
x = E/uiq [see Eq. (|18l) in the main text] and the imaginary part 7 = Gtn(ro)/27r. The latter appears due to the 
coupling to the tunneling tip. The terms ±.i-k/2 in Eqs. (|4"T)) . (|4"5jl appear due to the topological term Sg[Q], see Eq. 
© . With the above expressions for DQ and the action we can calculate the total partition function of the theory, 

J DQ e~ s[Q] = J tanh | dn dij d( exp 

This is exactly what one expects for the partition function of a supersymmetric theory. Note that the topologically 
nontrivial sector does not contribute to this result since the action 5*2 contains no Grassmann variables. 

In order to calculate the density of states, we need the pre-exponential factor from Eq. (1131) . In our parameterization, 
it has the following form in the two sectors: 

{cosh 2 -, Of = 0i, 

K K ( 50 ) 

sinh 2 ^-ir)tcosh 2 ~, 6 f = 2 . 



2i-kx I smh — — i?7<, cosh — 



VK 

~2 



= 1. 



(49) 



p = v cos Q\ Re / tanh — dn dr\ <i£< cosh — exp 



^sinh 2 — — irjC, cosh 2 — ) exp 



2iirx cosh — — — 



= v cos 6*i Re 



dn cosh — sinh — 
2 2 



2inx cosh 2 — exp (^2inx sinh 2 



v cos 0i 1 — Re 



1 + e 2 



2iitx 



v cos u\ 





17T 


2) + 










" ~2 




f 






2) + 


exp 



exp I 2iirx cosh 



n(x 2 + 7 2 ) 



sin(27ra;) 

2~KX 



(51) 



In the last expression we have used the condition 7 <C 1. The result (|STj) coincides with Eqs. (|T5|) - (fTT)) of the main 
text. In the tunneling limit 7^0 the lorentzian term yields a delta function. This implies that the Majorana state 
is robust in a closed system and disorder does not smear it. 

Note, that for an alternative model of the symmetry class D, when the topological term is absent, the density of 
states does not contain the Majorana delta-peak. Two topological sectors contribute with the same sign yielding [l5[ 



Pd = v cos Hi 



1 + 



sin(27rx) 



2ttx 



(52) 



High energies Exh <B< A 

In the previous section, we have calculated the density of states at low energies. The result (I5T1) exhibits a delta- 
peak at zero energy and oscillations which decay at the characteristic energy scale ljq. These oscillations appear due 
to repulsion between low-lying levels with energies E and —E. The scale ujq is the global mean level spacing inside 
the hole. At higher energies E >• w , the effect of level repulsion can be neglected and the density of states is given 
by the mean-field expression 

p(r,E) = i/Recos0(r,£), (53) 

where is a solution of the Usadel equation (fT0|) with k = 0. The result (|53[) can be obtained from the sigma model 
identity ([T3|) by integrating over Q around the unique global minimum of the action given by the Usadel equation. 
Supersymmetry ensures that the integration with respect to small fluctuations around the minimum yields unity and 
hence Eq. (|5"3")) follows. 

The Usadel equation is a complicated non-linear equation that cannot be solved analytically for arbitrary E. At 
relatively low energies luq <C E <C -Etii, the solution can be found perturbatively in E. While linear correction in E 
to is purely imaginary and does not change the density of states, the second order calculation up to terms ~ E 2 is 
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Figure 3: Spatial dependence of the density of states for different values of E. At lowest energies, the result (|16[) of the main 
text is reproduced, while at energies well above Ej-h the density of states is depleted only close to the hole edge within the strip 
in agreement with Eq. ()55[) . Remarkably, the density of states is first increased above its normal value before 
dropping to zero at r — R. 



necessary to obtain the observable result. The spatial profile of this small energy correction does not factorize as in 
the main term, see Eq. (I15[) in the main text. Second order perturbation theory amounts to solving linear differential 
equations, obtained by linearizing Eq. (|10p . and leads to very complicated functions. We will not discuss any further 
details of this approach here. 

As the energy is increased further and exceeds the Thouless energy, i?Th « £ « A, the Usadel equation admits 
another approximate solution. The angle 9 is small everywhere except a narrow vicinity of the hole boundary. In this 
vicinity we can neglect the curvature of the superconductor edge and write an approximate one-dimensional Usadel 
equation 



d 2 9 

D —z + 2iE sin 9 = 0. 
or z 

With the boundary conditions 9(R) = tt/2 and 9(r <C R) — > 0, solution to this equation reads 



(54) 



9 = 4 arctan 



( V2 - 1) exp 



-2i- 



E 



(- r 



(55) 



The density of states is given by Eq. (|53p with the above result for 9. We see that 9 decreases from tt/2 down to 
exponentially small values in a narrow strip of the width ~ R^J E^h/E near the hole boundary. The density of states 
takes its normal value v everywhere inside the hole except for this narrow strip where it is depleted down to zero. 
Spatial integration yields global density of states 



N(E > E Th ) = ttvR 2 



1 - (2 - V2 



■ETh 



(56) 



This is the result (|20l) of the main text. 

We have also solved the Usadel equation numerically for the whole range of energies below A. The spatial profile of 
the density of states is shown in Fig.[3]and the energy dependence of global DOS is depicted in Fig.f2]in the main text. 
It perfectly matches both the low-energy limit (up to delta peak and oscillatory term) of Eq. (fT5l) and high-energy 
asymptotics ([56| . 



Tunneling current and noise 



Tunneling current and noise can be found using the Landauer formalism instead of the sigma model. Electrons 
incident from the metallic probe are either reflected normally or experience Andreev reflection. Amplitudes of these 
processes can be found from the microscopic Hamiltonian of the system and used to calculate the current and higher 
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cumulants, such as noise. The current I(V) found this way coincides with the result derived from the sigma model 
and presented in the main text, Eq. (|2ip . The noise, Eq. (|23[) . was found by means of the Landauer approach since the 
sigma-model consideration of this problem is rather cumbersome. A detailed derivation of the Landauer formalism 
for our problem and the calculation of current and noise will be presented elsewhere [1(1] . 
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